active_sites = [1 2 4 7 8];
passive_sites = [3 5 6 9];
control_sites = 10;

for site = 1:9
    for layer = 1:4
        for time_point = 1:4
            cur_act = proj_meta(site).rd(layer,time_point).act;
            cur_fb_bi = zeros(size(cur_act(1,:)));
            cur_pb_bi = zeros(size(cur_act(1,:)));

            % generate binary activity during feedback.
            velM_fb_bi = proj_meta(site).rd(layer,time_point).velM_smoothed(1:7500);
            velM_fb_bi = -velM_fb_bi; % make the velocity of forward-running positive.
            velM_fb_bi(velM_fb_bi > 0.001) = 1;
            velM_fb_bi(velM_fb_bi <= 0.001) = 0;
            
            % generate binary activity during playback.
            velM_pb_bi = proj_meta(site).rd(layer,time_point).velM_smoothed(7501:15000);
            velM_pb_bi = -velM_pb_bi; % make the velocity of forward-running positive.
            velM_pb_bi(velM_pb_bi > 0.001) = 1;
            velM_pb_bi(velM_pb_bi <= 0.001) = 0;
            
            % generate the frames of ps onset.
            ps_ons = find(diff(proj_meta(site).rd(layer,time_point).ps_id(1:7500))>1);
            % remove early & late ps.
            ps_ons(ps_ons > 7425) = [];
            ps_ons(ps_ons < 6) = [];
            
            % generate the frames of playback halt onset.
            pb_ons = find(diff(proj_meta(site).rd(layer,time_point).ps_id(7501:15000))>1);
            % remove early & late pb.
            pb_ons(pb_ons > 7425) = [];
            pb_ons(pb_ons < 8) = [];
            
            rps_nbr = 0; % number of real ps
            rps_id = [];
            for ps_id = 1:length(ps_ons)
                if sum(velM_fb_bi(ps_ons(ps_id)-5:ps_ons(ps_id)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = ps_id;
                    all_rps_nbr(site,time_point) = rps_nbr;
                end
            end
            
            rpb_nbr = 0; % number of real pb halt
            rpb_id = [];
            for pb_id = rps_id
                if sum(velM_pb_bi(pb_ons(pb_id)-7:pb_ons(pb_id)+22)) == 0
                    rpb_nbr = rpb_nbr+1;
                    rpb_id(rpb_nbr) = pb_id;
                    all_rpb_nbr(site,time_point) = rpb_nbr;
                end
            end

            rps_resp_cnt = 0;
            rps_resp = [];
            for rps_act_nbr = 1:length(rps_id)
                rps_ons = ps_ons(rps_id(rps_act_nbr));
                cur_fb_bi(rps_ons+31:rps_ons+45) = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
            end
            
            rpb_resp_cnt = 0;
            rpb_resp = [];
            for rpb_act_nbr = 1:length(rpb_id)
                rpb_ons = pb_ons(rpb_id(rpb_act_nbr));
                cur_pb_bi(rpb_ons+7531:rpb_ons+7545) = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
            end
            
            mm_id = 0;
            nm_id = 0;
            ph_id = 0;
            np_id = 0;
            mm = [];
            nm = [];
            ph = [];
            np = [];
            for frame_id = 1:length(cur_act)
                if cur_fb_bi(frame_id) == 1
                    mm_id = mm_id+1;
                    mm(mm_id) = frame_id;
                end
                if cur_fb_bi(frame_id) == 0
                    nm_id = nm_id+1;
                    nm(nm_id) = frame_id;
                end
            end
            for frame_id = 1:length(cur_act)
                if cur_pb_bi(frame_id) == 1
                    ph_id = ph_id+1;
                    ph(ph_id) = frame_id;
                end
                if cur_pb_bi(frame_id) == 0
                    np_id = np_id+1;
                    np(np_id) = frame_id;
                end
            end
            
            mm_act = [];
            nm_act = [];
            ph_act = [];
            np_act = [];
            for cell_id = 1:size(cur_act,1)
                mm_act(cell_id,:) = cur_act(cell_id,mm);
                nm_act(cell_id,:) = cur_act(cell_id,nm);
                ph_act(cell_id,:) = cur_act(cell_id,ph);
                np_act(cell_id,:) = cur_act(cell_id,np);
            end
            yas1{site}(layer,time_point) = struct('mm',mm_act,'nm',nm_act,'ph',ph_act,'np',np_act);
        end
        cell_cnt{site,layer} = size(proj_meta(site).rd(layer,1).act,1);
    end
end

